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In computer simulations, quantum delocalization of atomic nuclei can be modeled making use of 
the Path Integral (PI) formulation of quantum statistical mechanics. This approach, however, comes 
with a large computational cost. By restricting the PI modeling to a small region of space, this cost 
can be significantly reduced. In the present work we derive a Hamiltonian formulation for a bottom- 
up, theoretically solid simulation protocol that allows molecules to change their resolution from 
quantum-mechanical to classical and vice versa on the fly, while freely diffusing across the system. 

This approach renders possible simulations of quantum systems at constant chemical potential. 

The validity of the proposed scheme is demonstrated by means of simulations of low temperature 
parahydrogen. Potential future applications include simulations of biomolecules, membranes, and 
interfaces. 

PACS numbers: 05.10.-a,82.20.Wt,05.30.-d,61.20.Ja 

INTRODUCTION 


Nuclear quantum delocalization plays a crucial role in 
low temperature systems, e.g. helium or hydrogen [1-6], 
which can undergo a superfluid transition, and it affects 
in nontrivial ways a large variety of systems and pro¬ 
cesses at more standard thermodynamic conditions. It is 
the case, for example, for proton transfer in biomolecules 
and membranes and in DNA oxidation [7-14], the ther¬ 
modynamics of ice [15], the structure of water adlayers 
on catalysts [16, 17], and the structure and dynamics of 
bulk water at room temperature [18-22]. In order to 
account for these effects in computer simulations, one 
can make use of Feynman’s Path Integral (PI) formula¬ 
tion of quantum statistical mechanics [2, 23-25], which 
enables the accurate description of nuclear delocaliza¬ 
tion by means of Monte Carlo (MC) or Molecular Dy¬ 
namics (MD) simulations [2, 24, 25]. This possibility, 
though, comes at the expense of an increased compu¬ 
tational cost. A strategy to overcome this limitation is 
to restrict the PI description of the atoms to a (small) 
region of space, where their quantum nature has to be 
explicitly accounted for and to model the other atoms as 
classical particles interacting via an appropriately chosen 
effective potential. Molecules diffusing across the bound¬ 
ary separating these two regions must change “on the 
fly” their representation from classical or quantum and 
vice versa. This approach, which is obviously viable only 
for sufficiently short De Broglie wavelength, is beneficial 
especially for applications where one has, at the same 
time, a small region which requires a PI description in a 
much larger simulation box. Textbook examples of such 
systems are given by liquid-solid or liquid-liquid inter¬ 


faces [26-28] and in protein simulations [29-31] in which 
a chemically accurate model of the active site can be con¬ 
currently employed with a coarser description of the rest 
of the molecule. The simplified model in the classical re¬ 
gion can also allow one to change on the fly the number 
of molecules in the system [32], thereby implementing 
a grand canonical PI approach. More generally, an ap¬ 
proach in which a classical and a PI model of the system 
are concurrently used in the same setup would allow a 
substantial computational gain. In turn, this enables the 
simulation of significantly longer length scales and sam¬ 
pling times compared to fully quantum PI simulations. 

A first step in this direction was taken in the frame¬ 
work of the Adaptive Resolution Simulation (AdResS) 
scheme [33-35] by merging quantum and classical effec¬ 
tive forces [36-38]. This work demonstrated the possi¬ 
bility of investigating the properties of a system of light 
atoms or molecules by explicitly considering their quan¬ 
tum nature only locally, without disrupting the overall 
thermodynamic balance between the quantum and the 
classical regions. However, the AdResS scheme is in¬ 
trinsically based on the interpolation of forces, and does 
not admit a Hamiltonian formulation [39]; therefore the 
quantum-classical coupling was introduced ad hoc after 
the quantization of the system and the introduction of 
fictitious momentum coordinates. This approach is thus 
incompatible with a proper PI quantization. 

In this paper, we provide a theoretically solid quantum- 
classical coupling protocol, based on a global Hamilto¬ 
nian. We simulate a system of atoms or molecules that 
exhibit quantum behavior only in a restricted region of 
space and behave as purely classical particles everywhere 
else. Additionally, we allow molecules to freely diffuse 
across the simulation domain and switch the nature of 
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FIG. 1. Illustration of the simulation setup for the quantum- 
classical simulations. Red (resp. blue) color corresponds to a 
larger (resp. smaller) radius of gyration. The smooth transi¬ 
tion from extended to collapsed molecules demonstrates the 
transition from quantum mechanical to classical behavior. 
The particles freely move between the regions and change 
their behavior accordingly. 

their interactions according to their position in space. 

ADAPTIVE QUANTUM-CLASSICAL COUPLING 

In order to model each subsystem with its appropriate 
interactions, we make use of the Hamiltonian AdResS 
(H-AdResS) method [40-42]. This scheme was devel¬ 
oped to perform adaptive resolution MD/MC simulations 
based on a global Hamiltonian, which makes it the ap¬ 
propriate framework for our work. A typical H-AdResS 
system is partitioned in two regions connected via a hy¬ 
brid buffer region. The resolution of a molecule depends 
on the value of a representative coordinate R (usually 
chosen to be the center of mass), and is parametrized by 
a continuous function A(R) smoothly switching from 0 
to 1 in the hybrid region. The total potential energy of 
each molecule is obtained by interpolating between the 
two resolutions. The H-AdResS Hamiltonian If of a sys¬ 
tem of point-like particles reads 

N 

R = A + E + (! - X «)v° - AH(X a )] (1) 

OL= 1 

where /C is the kinetic energy, a indexes the N particles, 
and A a = A(R Q ). The single-particle potentials V^ es 
(with Res = 0,1) are the sums of all intermolecular po¬ 
tentials acting on particle a , properly normalized so that 


double counting is avoided [40, 41]. In the following we 
make no assumption about the specific form of these in¬ 
teractions. The term AH, referred to as the Free Energy 
Compensation (FEC) [40, 41], is an external field acting 
only in the hybrid region to neutralize the density im¬ 
balance that naturally occurs when different models of 
the same system are coupled together. Its calculation is 
described in the Validation section. 

The employment of the H AdResS Hamiltonian in the 
PI formalism is straightforward. Specifically, the ring 
polymer potential energy obtained from the PI quantiza¬ 
tion of the Hamiltonian in Eq. (1), assuming Boltzmann 
statistics, is given by 

N P ( 2 

vp = E i 2 

a=l fc=1 (2) 

+ p [ X a,kV^ k + (1 — A a,k)V® k - AIf(A Qj fe)] j 

for N interacting particles in 3 dimensions, where up = 
VP/Ph, ft = 1 fksT, T is the temperature, kp is Boltz¬ 
mann’s constant, and h is Planck’s constant. The index 
k labels the P “copies” of the original system, which, 
after quantization consists of N ring-polymers, each con¬ 
taining P beads, X a , k = X(r a , k ), and V^ k s is the total 
interaction of type Res = 0,1 on replica k of particle a. 

Eq. (2) describes a system of quantum particles, rep¬ 
resented by ring polymers whose interactions change in 
space. Nevertheless, their quantum behavior, dictated 
by the strength of the springs connecting the beads of 
each ring, is the same everywhere. At this stage we need 
a strategy to switch between the quantum and classical 
descriptions of the particles. This can be achieved by 
modifying the mass of the atoms, as larger masses corre¬ 
spond to stronger springs of elastic constant mup ; a large 
mass causes the ring polymers collapse, and the particles 
approach their classical limit. We thus define 

m -A fi( A) = Am + (1 — A )M (3) 

where /z(A) smoothly switches from a mass p,(0) = M 
to a mass /i( 1) = rn <C M. For fi = m, which is set 
to be the real, physical mass, the particles are light and 
the quantum zero-point motion becomes important. In 
contrast, the mass M should be large enough to give the 
particles a classical character. 

We now proceed with the quantization of a system of 
particles with position-dependent masses. As a starting 
point we consider the Hamiltonian operator for a free 
particle of mass /x( x) in one dimension (the procedure 
generalizes trivially to many-particle systems in three di¬ 
mensions). The Hamiltonian must be represented as a 
Hermitian operator, which we can obtain by writing it in 
the following form: 

Td=^P^ 1 {x)p (4) 
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where p- x {x) is the inverse mass operator, and p is the 
momentum operator. Using this Hamiltonian, we seek to 
formulate the partition function Q = Tr[exp{— pH}] as a 
path integral. Introducing the usual set of P resolutions 
of the identity operator, we can write the trace as 


Q = lim 

P—> OO 


dx i • ••dxp 


p 

\ (x k |exp 
k =1 




\x k + 1 ) 


X P + 1=X 1 


( 5 ) 


Note that in the free particle case, it is not necessary 
to make use of the limit P — > oo, required when applying 
Trotter’s theorem. However, the latter is generally neces¬ 
sary in presence of a potential V(x), hence we introduce 
this limit at this stage without any loss of generality. 

To derive the matrix elements in Eq. (5), we introduce 
the momentum identity resolution: 


{x k \ exp ( pn 1 (x)p ) \x k+l ) 


J dp(x k \p) (p|exp \xk+i) 


( 6 ) 


ing eigenvalues. Substituting Eq. (9) into Eq. (6) gives 


(xk | exp (~jpPP 1 {x)p S j \x k+1 ) 

/•OO 

= / dp(x k \p) {p\x k+ i} x 

P 


x exp 


= ( p(x k+1 )P \ 

\ 27 xj3h 2 ) 


2 P 

1 

2 

exp 


P p (x k+ i )+ihp 


dp 


-l 


dx 


%k + l / 


/3p(x k +i)P 

2 m 2 


(x k - x k+ i)+ 


/3H 2 dp 1 


2 P dx 



2 1 

%k + l- 

j 


( 10 ) 


where the last equality has been obtained by introducing 
the matrix elements (x\p) = exp(ipx/ h) / y/2nh and per¬ 
forming the momentum integration by completing the 
square. 

From Eq. (10), we see that the inverse mass derivative 
term can be neglected if the following condition holds: 





2-PAxfc^^.1 
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( 11 ) 


Given that the limit P — > oo is ultimately taken, we 
can work with an infinitesimal version of the exponential 
operators by expanding the exponential to first order. 
Thus, we obtain 


where we defined A.x k , k +i = \x k — x k + 1 |. Since, 


(dp- 1 ' 

\ 


1 

( dp\ 

\ dx , 

Xk + 1 


P 2 (x k+1 ) 

\ dx Jx k+1 


(p\ exp ( --j^ppp \x)p)\x) 
{P\ ( 1- jpPP~\x)P ) |x> 


( 7 ) 


Now, we introduce the commutator [p 1 (x),p] and write 
p~ 1 {x)p = pp~ 1 (x) + [// _1 (x),p] 


= pp 1 ( x)+ih 


dp 


-l 


dx 


( 8 ) 


Substituting Eq. (8) into Eq. (7) yields 


, | (. P ~ 2 -i/m < h3 J l P 1 \ | \ 
/ \ \ (^ Pp 2 -u \ ih P d P~ 1 \ 

= <pk)(i-jpc M-2 pP-^) 


( p\x ) exp 


P ( 2 - 1 ,_\ , 


-l \ -i 


2 P \ P "" W +ihp -dx 


( 9 ) 


where the operators are now replaced by the correspond- 


the condition becomes 





2^x kk ^i 

A l(xk+i) 


p(x k+ i) 


(13) 


using the definition of position-dependent De Broglie 
wavelength A M (x) = \J/3h 2 /(P p(x)) . Since (Ax) = 

\/(FEf=iH+i) = A MV / ( p - 1 )/ p ~ for a free 

ring of constant mass p and typical values of P, we can 
approximate Ax/^fc+i ss A At (xfe+i) and write 


dp(x) 2p(x) 

dx A M (x) 


(14) 


for an arbitrary position x. The inequality in Eq. (14) 
must be satisfied everywhere in the system. In the clas¬ 
sical and quantum regions this is trivially the case, as 
the resolution function A(x) is flat there. This condition 
means that the interpolation within the hybrid region 
needs to be sufficiently smooth in order to neglect the 
term containing the mass gradient in Eq. (10). This 
can always be achieved by utilizing a sufficiently large 
coupling region. In that sense, the criterion can be inter¬ 
preted as a lower bound on the width of the hybrid region. 
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Furthermore, it also holds in the presence of typical po¬ 
tentials, since these typically do not dramatically change 
the radius of gyration and the intrabead distances of the 
polymer rings compared to free rings. Additionally, al¬ 
though the derivation was carried out for one dimensional 
systems, the derivation generalizes trivially to higher di¬ 
mensions. This criterion is also correct in three spatial 
dimensions, as the mass change only happens along one 
of these dimensions, and therefore, only the bead-bead 
distances projected onto this direction matter. 

Concluding, if the inequality in Eq. (14) is fulfilled, 
then introducing the H AdResS potential energy, we ob¬ 
tain the following partition function for N interacting 
Boltzmann particles in three dimensions: 


P—yoo 


0 -- ta nn 


.k—1 a=1 
» P N 


nn dr a ,k exp {~PVp} 


k—1 a—1 


(15) 


pronounced quantum mechanical character [4-6]. We 
hence study it as an extreme case, well-suited to test 
the proposed algorithm. 


System setup 

We consider a system composed of 4964 hydrogen 
molecules in a slab of dimensions 24.000 nm x 3.123 nm x 
3.123 nm (molecular density 28.4cm 3 /mol) with periodic 
boundary conditions in all directions. The width of 
the low-mass quantum region is set to c?qm = 6.0 nm 
and the thickness of each hybrid transition region is 
dnY = 5.0 nm. In order to assign to a bead its position- 
dependent resolution A, its distance from the boundary 
between the quantum and the hybrid region is computed, 
i.e. \x a ^k\ — dQM/2, where x a> k denotes the X coordinate 
of the bead in a coordinate system with its origin at the 
center of the simulation box. This quantity is then em¬ 
ployed in the resolution function A(x), which is given by 


with 


N 


*? = E E 


k= 1 Q! = l 


f-Lot.k | |2 ^ i fJ'cx.k 

—2~ Ir^-r^l -^ l0 § — 


+ p + (1 — ^a,k)V^ k — Aff(A a ,fc)] | 


(16) 


and n at i = n(v a ^). In Eq. (16) the position-dependent 
normalization prefactor has been explicitly introduced in 
the potential Vp as a logarithmic function of the bead 
masses, so that it can be treated as a conventional energy 
term and fully removed from the Hamiltonian by means 
of the FEC function AH in Eq. (1), in a manner similar 
to that done in Ref. [43]. The light mass m has been used 
as the reference mass scale. A different choice would not 
affect the final result of the calculations. Using the mass 
m as a reference, however, the normalization prefactor 
corresponds to the one known for Pis with constant mass 
m [23, 25]. The ring polymers described by the energy 
function Vp (Eq. (16)) are expanded in the region where 
the mass is small and collapse to nearly classical point- 
like particles in the large-mass region. 


VALIDATION 

To validate the proposed quantum-to-classical cou¬ 
pling scheme, adaptive Path Integral MC simulations of 
liquid parahydrogen at 20 K with P = 16 are performed. 
Other test cases might be considered, e.g. water at room 
temperature, but, in spite of the important role played 
by nuclear quantum effects in this example [18-22], the 
hydrogen atoms feature a relatively small delocalization. 
Ultracold hydrogen, on the other hand, exhibits a more 


{ 1 : x < 0 

cos2 (f dTx) '■ 0 < x < day (17) 
0 : x > c?hy 


The mass m is set to the molecular hydrogen mass 
tor 2 = 2.001 au. In the classical region the increased 
mass is chosen as M = 100toh 2 - In the quantum (QM) 
region we employ the Silvera-Goldman potential [44, 45] 
with a cutoff at 0.9 nm for the intermolecular interaction 
potential V 1 , while in the classical (CL) region we make 
use of a shifted, purely repulsive Weeks Chandler Ander¬ 
sen (WCA) potential [46]: 


U° 



if r < R c 
if r > R c 


(18) 


where r = \r a ^k — »"/ 3 ,fc|, (a ^ /3) denotes the distance 
between beads of the same imaginary time slice k in 
different molecules a and /?. Furthermore, we choose 
e = l.OkJ/mol, a = 0.14nm, and r$ = 0.15nm. The 
cutoff is given by R c = 2^ 6 a + Tq. The two potentials 
are graphically presented in Fig. 2. This potential is 
not to be interpreted as a classical model for low tem¬ 
perature parahydrogen; rather, it was parametrized only 
to approximately reproduce the hard-core radius of the 
reference quantum particles. Other choices, suitable to 
other simulation setups, are clearly possible. We pur¬ 
posely avoid fitting the classical potential to the struc¬ 
ture of the reference to demonstrate the generality of the 
protocol. 

The chosen set of parameters also satisfies Eq. (14). 
Finally, we stress that in the CL region the WCA in¬ 
teraction between ring polymers is computed only using 
the center of mass of the ring, thus gaining an effective 
reduction of the computational cost. This simplification 
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FIG. 2. Non-bonded intermolecular interaction potentials 
used in the adaptive quantum-classical simulations. The red 
curve is the Silvera Goldman potential, which is employed 
in the low-mass quantum region. The blue curve shows the 
shifted WCA potential, which is used in the high-mass clas¬ 
sical region. 

is allowed by the essentially point-like structure of the 
rings in the CL regions, as can be seen from the radius of 
gyration profile (Fig. 3). The number of computations 
per pair of molecule is reduced from P = 16 to one. 

To modulate the thermodynamic imbalance between 
the classical high-mass and the low-mass quantum sub¬ 
systems a Free Energy Compensation (FEC) is ap¬ 
plied [40, 41]. To compute the AF/kti(A) compensa¬ 
tion a Kirkwood Thermodynamic Integration [47] of a 
smaller system of 360 molecules in a box with dimen¬ 
sions 2.570 nm x 2.570 nm x 2.570 nm is performed. 

In order to remove also the remaining fluctuations in 
the obtained density profile after applying the Kirkwood- 
based FEC, an iterative approach similar to the one pre¬ 
sented in [35] is employed. The normalized density profile 
Pwy(x) hi the hybrid region is transformed into a func¬ 
tion of the resolution, phy(A), and then converted into a 
correction energy of the form 

Atf(A) = —iln{p H Y(A)} (19) 

The latter quantity is then applied as part of the FEC 
AH(X) in addition to the term AUkti(A) obtained from 
Kirkwood thermodynamic integration. This is done in an 
iterative fashion, until a sufficiently flat density profile is 
obtained. The protocol for the FEC then reads 

AH i+1 (X) = AF(A) - iln{pj IY (A)} (20) 

with AH°(X) = AF/kti(A) and PhyM = Phy(^)i 
where Ppy^A) corresponds to the initial hybrid region 


density profile obtained from simulations in which only 
the Kirkwood-based FEC term is applied. The protocol 
converges by construction when a flat density profile is 
achieved. 

Monte Carlo sampling 

To sample the system’s phase space we employ a stan¬ 
dard Metropolis Monte Carlo algorithm [25]. For the 
Kirkwood TI of the small system we run 16 simulations 
with 10 5 sweeps each. The A parameter increases linearly 
every sweep by 10 -5 . The results are averaged after the 
simulations. Employing the Kirkwood TI FEC term thus 
obtained we then run 5 iterations of simulations with ap¬ 
plying the protocol set out above to refine the density 
profile. Each iteration consists of 32 parallel simulations 
with each of these running 5-10 3 equilibration sweeps and 
another 5 • 10 3 sweeps during which the density profile is 
measured. Also here, after each iteration the results are 
averaged. Having reached a sufficiently smooth density 
profile, we then utilize the FEC from the Kirkwood TI 
and the iterative protocol to perform the main produc¬ 
tion simulations. For these we perform 32 simulations 
in parallel, each running 4 • 10 3 sweeps. Afterwards, the 
results (i.e. the RDF’s, the density profiles as well as the 
radius of gyration profiles) are once again averaged over 
all simulations. 

Each sweep is constituted by N attempted Monte 
Carlo moves on randomly chosen molecules, with N being 
the total number of molecules (N = 4964 in the produc¬ 
tion run simulations). Three different kinds of moves are 
randomly performed: 

Whole molecule displacements: The chosen 
molecule is displaced as a whole by moving its center of 
mass. The direction is chosen randomly from a uniform 
spherical distribution and the distance is drawn from a 
Gaussian distribution with zero mean and width crcoM- 

Molecule rotations: The chosen molecule is rotated 
as a whole around a randomly oriented axis passing 
through its center of mass. The angle is chosen randomly 
from a Gaussian distribution with zero mean and width 

^rot- 

Individual Trotter-bead moves: An individual 
bead of the molecule is randomly chosen and displaced. 
The direction is chosen randomly from a uniform spher¬ 
ical distribution and the distance is drawn from a Gaus¬ 
sian distribution with zero mean and width crbead- 

The different values for the cq’s of all simulations are 
presented in Tab. I. In the adaptive resolution simula¬ 
tion, they are chosen such that they result in adequate 
acceptance ratios for the moves both in the classical high- 
mass as well as in the quantum low-mass region. When 
picking a molecule for a Monte Carlo move the proba¬ 
bilities for performing whole molecule displacements or 
molecule rotations are Y'l 3 each while the probability for 
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Trotter-bead moves was n /i3. This choice leads to a 
convenient balance between whole molecule motions and 
Trotter-bead fluctuations. 


Simulation 

& CoM 

&rot 

^"bead 

Kirkwood TI 

0.1 nm 

0.5 rad 

0.03 nm 

Adaptive Simulation 

0.1 nm 

0.5 rad 

0.03 nm 

QM Reference 

0.1 nm 

0.5 rad 

0.07 nm 

Classical Reference 

0.1 nm 

- 

- 


TABLE I. Widths of the Gaussian distributions employed to 
draw the random displacements and rotations from for the 
Kirkwood Thermodynamic Integration, for the reference sim¬ 
ulations as well as for the adaptive quantum-classical simula¬ 
tions. 


Reference simulations 

To be able to evaluate the results of the adaptive 
quantum-classical simulations, we perform full-quantum 
as well as full-classical reference simulations of liquid 
parahydrogen for comparison. The systems are composed 
of 828 molecules in a box with dimensions 4.003 nm x 
3.123 nmx 3.123 nm. These parameters result in the same 
density as in the adaptive simulations. Likewise, the tem¬ 
perature is set to T = 20 K and the Silvera-Goldman 
potential is employed. For the full-quantum simulations 
we choose P = 16 as in the adaptive simulations while 
the classical simulations are performed with P = 1. In 
both cases, 16 simulations are run in parallel, each one 
for 2 • 10 4 sweeps. Afterwards the results are averaged. 
The values used for the cq’s in the reference simulations 
are presented in Tab. I. In the classical simulations, all 
moves are as described above with the obvious excep¬ 
tion of bead and rotating moves, which do not exist for 
classical particles. 

Results 

A snapshot of the dual-resolution simulation is pre¬ 
sented in Fig. 1: the gradual change in size of the ring 
polymers indicates the transition from the classical to 
quantum mechanical regions and vice versa. Results are 
reported in Fig. 3. The radius of gyration of the ring 
polymers in the quantum region (QM) perfectly repro¬ 
duces the one of a corresponding fully quantum simula¬ 
tion. In the CL region the radius of gyration drops by 
ss 90%, indicating the classical character of the molecules 
(see Fig. 1). By means of the FEC a nearly flat density 
profile was obtained in the quantum region. 

A quantitative measure of the fluid structure is pro¬ 
vided by the radial distribution function (RDF): the lat¬ 


ter is computed only in the inner part of the QM region 
(shaded region in Fig. 3). In spite of the remarkable 
differences between the quantum fluid and the classical 
model, the RDF measured in the QM region matches very 
well the one obtained in the completely quantum refer¬ 
ence simulations. These results show that in the QM re¬ 
gion the quantum-to-classical coupling scheme correctly 
reproduces the structure of the quantum mechanical sys¬ 
tem. 
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FIG. 3. Top: profiles of the radius of gyration r g and the 
normalized density p. The blue reference corresponds to the 
radius of gyration of molecules in a corresponding completely 
quantum reference simulation. The shaded area marks the 
region used for the calculation of the RDF. Bottom: RDFs 
of the quantum-to-classical simulation calculated in the in¬ 
ner part of the quantum region (red), of a corresponding 
completely quantum reference simulation (blue), and of a 
full-classical system of particles interacting via the Silvera- 
Goldman potential (green). 


Speedup over full-quantum simulations 

As mentioned earlier, in the proposed quantum-to- 
classical coupling scheme, interactions in the classical re¬ 
gion do not need to be calculated P times, with P being 
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the Trotter number, but because of the collapse of the 
polymer rings they are computed only once between the 
centers of mass of the (quasi point-like) rings. Addition¬ 
ally, a numerically simpler potential with a shorter cutoff 
is used in the classical region. Therefore, simulations be¬ 
come computationally more efficient, as we demonstrate 
hereafter. 

Since in practice the method is most beneficial for sys¬ 
tems in which the classical region is much bigger than 
the quantum region, we will consider such a situation. 
We perform 4 sets of simulations for different box sizes 
with each set consisting of a full-quantum, a full-classical 
and adaptive simulation in which the quantum region has 
a width of 2.0nrn and the adjacent hybrid regions each 
have widths of 1.0 nm. The total box sizes as well as the 
corresponding molecule numbers for the simulations are 
presented in Tab. II. The temperature and density are 
the same as before. In all cases, the classical regions 
are significantly larger than the quantum ones. Note 
that, here, “classical simulation” denotes a simulation of 
a WCA-liquid of collapsed polymer rings, exactly as in 
the classical region of the adaptive simulations, and not 
of a classical, i.e. with P = 1, liquid of parahydrogen. 
All simulations are run for 400 sweeps and in the case of 
the adaptive simulation, a FEC is applied. Furthermore, 
the set of Monte Carlo moves is chosen for all of them to 
be that previously used for the adaptive simulations. 


Number of molecules 

Lx 

Ly 

L z 

4964 

24.0 nm 

3.123 nm 

3.123 nm 

6619 

32.0 nm 

3.123 nm 

3.123 nm 

8273 

40.0 nm 

3.123 nm 

3.123 nm 

9928 

48.0 nm 

3.123 nm 

3.123 nm 


TABLE II. Number of molecules and box geometries for the 
different sets of simulations for the calculation of the compu¬ 
tational gain of adaptive quantum-classical over full-quantum 
simulations. 

In general, the overall speedup of the simulations 
strongly depends on the details of the implementation 
of the algorithm and is therefore platform dependent. 
For example, if there is a high overhead in the code, 
the overall computational gain by more efficient poten¬ 
tial energy calculations will be small. If the program 
spends most of its time with these calculations, a sig¬ 
nificant improvement is possible. Hence, in order to ob¬ 
tain platform independent results, we only measure the 
time our code spends with potential energy calculations. 
These times are plotted in Fig. 4. Additionally, the cor¬ 
responding speedups, defined as T quantum /T ada pti V e with 
^adaptive (^quantum) being the time spent for the energy 
calculations in the adaptive (quantum) simulations, are 
presented. 

It can be seen that the adaptive simulations are sig- 
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FIG. 4. Main figure: times required for the potential energy 
calculations in quantum, classical and adaptive simulations 
for systems of four different box lengths L x . Inset: speedup 
of the adaptive quantum-classical simulations. The lines are 
a guide to the eye. 

nificantly faster than their corresponding fully quantum 
counterparts. For the largest box, the energy calculations 
in the adaptive quantum-classical simulations are faster 
by a factor of ~ 4.5 than the full-quantum simulations. 
Furthermore, it is visible that the time required for the 
adaptive simulations stays nearly constant for the differ¬ 
ent box sizes. The reason for it is that the computational 
cost of the interactions between classical molecules in the 
larger simulations is negligible compared to the time re¬ 
quired for the computation of the potential energies in 
the quantum region. 

CONCLUSIONS 

In conclusion, we have derived a bottom-up 
Hamiltonian-based path integral formulation of a system 
of atoms or molecules whose quantum character depends 
on spatial location, and smoothly changes as the atoms 
or molecules diffuse. The formalism is derived with the 
aim of providing a new approach for treating quantum 
condensed-phase soft-matter problems at multiple levels 
of resolution, here, employing both quantum and “classi¬ 
cal” regions. Possible future applications are diverse and 
include, for example, adaptive quantum-classical simu¬ 
lations of interface systems, membranes, and proteins. 
The approach will also allow rigorous treatment of the 
quantum grand-canonical ensemble. Due to the reduced 
number of degrees of freedom in the classical subdomain, 
the protocol presented enables a computationally more 
efficient sampling of configurations compared to a fully 
quantum simulation. This, in turn, allows an extension 
of the accessible time- and length-scales. Furthermore, 












the proposed scheme can also be employed in more ad¬ 
vanced PI simulation techniques, such as Centroid Path 
Integral MD [25, 48, 49] or Ring Polymer MD [25, 50]. 
The development of such applications is the goal of a 
future study. 
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